      subroutine Solution
	    use global
        Implicit none
	    integer :: j,IPIV(nnt),INFO,i
	    real(8) :: M_L(nnt),Uj(nnt),Rj(nnt),DU_L(nnt),DU_H(nnt),M_n(nnt,nnt)

	    !!! Get the lump mass matrix.
	    M_L=sum(M,1)

	    do j=1,3
	      Uj=U(:,j)
		  Rj=Resid(:,j)
		  
		  !if (j==2) then
		  ! do i=1,nnt
		  !   write(*,'(3es15.4)') Uj(i),Rj(i),F(i,j)
		  ! enddo
		  !endif
		  
		  !!! Find the 1st order solution.
		  DU_L=(Rj+Cd*(MATMUL(M, Uj)-M_L*Uj))/M_L
		  !DU_L=(MATMUL(M, Uj)-M_L*Uj+Rj)/M_L
		  !if (j==3) then
		  ! do i=1,nnt
		  ! write(*,*) DU_L(i)
		  ! enddo
		  !endif
		  
		  !!! Find the 2nd order solution from [M]*[Rj]=[Uj].
		  !M_n=M  !! DGESV function will change the input mass matrix [M_n]. So it must be reinput.
		  !CALL DGESV( nnt, 1, M_n, nnt, IPIV, Rj, nnt, INFO )  
		  !!write(*,*) info
		  !IF (info==0)then
		  !  DU_H=Rj
		  !else
		  !  write(*,*) 'Error: can not find the 2nd order solution.'
		  !  exit
		  !endif

		  !!! Shock capture using FEM-FCT method.
		  !Call FEM_FCT
		  
		  !!! The solution increment from step 'n' to 'n+1'.
		  DU(:,j)=DU_L !+0.5D0*(DU_H-DU_L)
		  
	    enddo
	    
	    !!! Boundary condition
	    DU(1,2)=0.D0       ! The veloctiy of the fixed end of the chamber is zero. 
	    DU(nnt,2)=0.D0    ! So the solution increment is zero.
	    
	    !!! Update the solution vectors at step 'n+1'
	    U=U+DU
      end subroutine
